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Abstract 

We solve a problem of non-convex stochastic optimisation with help of simulated annealing of Levy 
flights of a variable stability index. The search of the ground state of an unknown potential is non-local 
due to big jumps of the Levy flights process. The convergence to the ground state is fast due to a 
polynomial decrease rate of the temperature. 

Keywords: Levy flights; simulated annealing; non-local search; heavy-tails; variable stability index; 
stable-like process; global optimisation. 



1 Introduction 

Let U be a potential function in M. d having several local minima and increasing fast at infinity. We look for 
a global minimum of U. Classical continuous-time simulated annealing (Boltzmann machine) (see |VG84[ 
IGidi IGH861 IGM93] ) consists in running a diffusion process 

Z 0z (t)=z-[ VU{Z 0z {u))du+ [ &{u)dW{u), 

U2 (1-1) 

where W is a standard Brownian motion, 9 > denotes the cooling rate and A > 1 parametrises the initial 
temperature, which equals y/6/ ln(A) at time t = 0. It is known that there is a critical value 9 such that 
the diffusion Z(t) converges in distribution to the global minimum of U if 9 > 9, and the convergence fails 
otherwise. Moreover, the critical value 8 is the logarithmic rate of the principal non-zero eigenvalue Ai(a) 
of a time homogeneous diffusion generator A a f = ^-A/ — (VU, /), i.e. 

9 = - lim cr 2 ln|Ai(cr)|. (1.2) 

The value of 9 can be calculated explicitly, if one knows the heights of potential barriers between different 
wells of U (see [Wen72al IWen72b] and [KM96] for precise results). Rigorous results on optimal cooling rate 
in simulated annealing algorithms can be found in |Haj88[ ICHS87) IHKS891 IHS90j . 

In order to accelerate the search, Szu and Hartley in |SH87j suggested the so-called fast simulated anneal- 
ing (Cauchy machine), which is a combination of a classical Metropolis algorithm introduced in |MRRT53] 
and a concept of non-local search due to the heavy-tail Cauchy visiting distribution. The authors claimed 
that in the Cauchy machine the temperature can be chosen decreasing as a power of time, namely a(t) ~ t , 
and applied the algorithm in image processing, see |Szu93j . 

Motivated by |SH87j , in our papers |Pav06b|, IPav06a| we considered a continuous-time counterpart of 
the process Z driven by Levy flights of stability index a £ (0,2), and temperature a(t) ~ t~ e , 9 > 0. We 
discovered that such a jump process never settles in the neighbourhood of a global minimum, but can be 
used to reveal a spatial structure of the potential U. The dynamics of Levy flights with constant small noise 
was studied in our previous papers [IP0 6a. IP06b], ITP06c] . 
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In the present paper we solve the problem of global optimisation with help of state- dependent Levy flights 
in a multi-well potential U. We show, that in certain annealing regimes, the global minimum is localised 
always, as in the classical Gaussian case. For simplicity, we restrict our theoretical argument in Sections ®- 
[5] to onc-dimcnsional potentials. However, it will be clear from the presentation, that the algorithm also 
works in a multi-dimensional setting. In our numerical example in Section \E\ we consider a two-dimensional 
potential with five wells. 

2 Results on the cooled down Levy flights 

In [Pav06b( lPav06a1 we considered a one-dimensional Levy flights process in an external potential U deter- 
mined by the equation 

= z - jf* U>{Z${u-)) du + jf* (2.1) 

We understand a Levy flights process as a symmetric stable Levy process with stability index a € (0, 2), 
whose marginal distributions have the Fourier transform 

Ee- i(Q) W = e-^M" = cxp (t [ [e**> 1 - iuyI D (y)] , (2.2) 

\ jr\{o} \y\ ) 

where Io(y) is the indicator function of the unit disk D = {y : \y\ < 1}, and c(a) = 2| cos (^) r(— a)\. We 
choose such a paramctrisation in order to have a simple form of a Levy measure v{dy) = \y\~ 1 ~ a dy. 

The measure v is also called a jump measure of the process Levy L^ a \ It controls the intensity and sizes 
of its jumps. Indeed, let AL^{t) = L^ a \t) - L<»(i-) denote the jump size of L" at time instance t > 0. 
Then the number of jumps on the time interval (0,i] with values in a set J C R is a Poissonian random 
variable with the mean tv(J) (which can be possibly zero or infinite). 

The process is a Markov process with a non-local generator 

Af(x)= [ \f(x + v)-f(x)-vf'(x)In(v)] T %;, (2.3) 

Jr\{o} \y\ + 

which is also referred to as a fractional Laplacian, A = A"/ 2 . 

We direct reader's attention to the books [Sat99[ |App04( IPro04j on a rigorous mathematical theory 
of Levy processes and stochastic differential equations. Physical results on the subject can be found in 
[MK041 ICGK+041 IBS02] . 

We assume that the potential U has n local minima mj and n—\ local maxima such that —00 = sq < 
mi < Si < ■ ■ ■ < m n < s n+ i = +00. The extrema are non-degenerate, i.e. U' '(mi) > and U"(si) < 0. 
Moreover we demand that |Z7'(a:)| > |x| 1+c as \x\ — * +00 for some positive c. 

In the small temperature limit, i.e. when A — > +00 or t — > +00, the process can be seen as a random 
perturbation of a deterministic dynamical system 

X° x (t)=x- f U\Xl{s))ds. (2.4) 
Jo 

We denote f2j = (s,_i, si), 1 < i < n, the domains of attraction of the stable points m,. 

The positive parameter 9 is called the cooling rate, and A > determines the initial temperature of the 
system which equals to A~ e at t = 0. 

Equation (|2.1[) describes a non-linear dynamics of a Levy particle, whose temperature is being decreased 
at a polynomial rate as t — > 00. 

In |Pav06bl[Pav06a1 . we discovered two cooling regimes — slow cooling 9 < 1/a and fast cooling 9 > 1/a 
— in which the transitions of a particle between the wells of U have different asymptotic properties. 

Let A > be a small number, and let Bi = {y : \y — mi\ < A} denote a A-neighbourhood of a local 
minimum m.;. Consider transition times 

T s % - hfi> > s : Z^(u) 6 U^fl,-} (2.5) 
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between different neighbourhoods Bi and the corresponding transition probabilities P SjZ (2^ Q '(T' i ' A ) G Bj), 
i 7^ j. Then for < 1/a and z G i?i we have 



Po, 2 (Z (Q) (T') G B 3 -) 



^l" 1 , A — > +00, 



«L a) [«f B) ]- 1 . Mi, 



where 



» 



» 



n 3 

dy 



\Sj-i - m, 



R\n 4 K-2/l 1+Q 



1 



1 



i ^ j, 
1 



a \ s»_i - mi 



(2.6) 
(2.7) 

(2.8) 



We have also shown that in the limit t — > +00, Z^,{t) has a distribution 

n 

7r^(dy)=J24 a) S mi (dy) 



i=l 



(a) 

1, Qu = ~<li 



(2.9) 

(a) 



It 



where the vector = (ir[ a \ . . . ,tt^) t , solves the equation Q T -n:^ = 0, Q = 

is clearly seen, that all 7r,[ Q ' ) > 0, and does not settle down near the global minimum of U. However, 

the values 7r l - Q ' ) , which can be estimated from the Monte Carlo simulations, reveal the spatial structure of U, 
e.g. the sizes of the domains fi^. 

If the cooling rate 6 is above the threshold 1/a, the Levy particle gets trapped in one of the wells 
and thus the convergence fails. Consider the first exit time from the i-th well 



Then, for z G Bi, 1 < i < n, 



and consequently, Eq, z S 1 = 00. 



Po, z {S l <oo)=0 



(2.10) 
(2.11) 



3 Levy flights with variable stability index (stable-like processes) 

In order to take into account the energy geometry of the potential, we have to make the Levy flights process 
depend on its current position. Thus instead of Levy flights L^ a ' defined in (|2.2p , we consider now the 
so-called stable-like process H = (H(t)) t >o, which is a Markov process defined by the non-local generator 

Bf(x)= [ [f(x + y)-f(x)-yf'(x)I D ( y )}—^--, (3.1) 
Jr\{0} |y|!+«W 

with a function a(x) taking values in the interval (0,2). Sometimes, the notation B = A Q (')/ 2 is used. 
The main difference between and H consists in a dependence of a stable-like jump measure v x [dy) = 
\y\~ l ~ a ( x *) on the spatial coordinate x. Thus, if H(to) = x$, the instant jump distribution of H at time to is 
governed by a stable measure v Xo (dy). 

Of course, the dynamics of H is completely determined by a variable stability index a(x). For example, 
if a(x) =Q(j£ (0, 2), then H is just a usual Levy flights process of index ctQ. From now on, we assume that 
a(x) takes values strictly between and 2, i.e. < a < a(x) < A < 2, to exclude degeneration of the jump 
measure. 

We are going to study the dynamics of a stochastic differential equation with the driving process H, 
namely 

Y 0ty (t)=y- fu'{Y^{u-))du+ ( t dH ^oAu-),u) ^ 

Jo JO \ A ' u ) 
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For a better understanding of the process Y, it is instructive to consider a discrete time analogue of (|3.2|) 
given by the recurrent formula 



y*=y^- u '^ h+ { xit%r fc - L (3 - 3) 

The discrete time dynamical system (|3.3j) is obtained from the Euler approximation of (|3.2p with the time 
step h, and can be used for simulations. (However, one should be careful when U' is not globally Lipschits.) 
The random input is determined by the random variables such that 

$(y) £ = fcVa(»)£(«(»))(l) (3.4) 

where L^ a ^\l) has a standard symmetric a(j/)-stable distribution with the Fourier transform 



4 One-well dynamics. Transitions and trapping 

The dynamics Y of the Levy flights with variable stability index in the force field U' is a result of an interplay 
of two independent effects. First, for small temperatures, i.e. when A — * +00 or t — ► +00, Y is close to the 
underlying deterministic trajectory X . Starting from any point of fit, it follows X° with the same initial 
point and with high probability reaches a small neighbourhood of m, in relatively short time. On the other 
hand, Y tries to deviate from X° making jumps controlled by the jump measure v x (dy). Finally, if Y is in 
the well f2j, it spends most of the time in a neighbourhood of m, and thus has jumps approximately governed 
by the stable jump measure v mi (dy) = \y\~ 1 ~ al - mi >dy. 

Thus, the exit time and the exit probability from the well Qi of the process Y are approximately the same 
as for the process Z( a ( m *)). This resemblance becomes exact if we consider a piece- wise constant stability 
index a(x), a(x) = Y^i=i a $-{ x £ < on < 2. With this choice of a(x), the process Y is just driven by 
the equation (|2.1[) until it exits the well. (We omit a discussion on the behaviour of the process in the small 
neighbourhoods of the saddle points.) 

Let us introduce the following transition and exit times for the process Y: 

t\ z = mf{u > s : Y SiZ (u) G Uj^Bj}, (4.1) 
<rl z = inf{u > s : Y s , z $ CU}. (4.2) 

It follows from (|2.6p and (|2.7p . that if y G Bi and a{rrii)6 < 1, then the following relations hold as A — >• +00: 



PoAY(T*) G Bj) q^\^ mi)) }-\ i + j. (4.4) 
On the other hand, if a{rrn)6 > 1, the Levy particle gets trapped in the well due to (|2.11|) . i.e. 

?oA<r i <°o) = \^+™, (4.5) 



and consequently, Eo^c 1 ' = 00. 

Relations ()4.3p . (|4.4p and (|4.5p are crucial for our analysis. 



5 Non-local random search and simulated annealing 

5.1 Getting trapped in an assigned well 

Wc demonstrate now how to drive the Levy particle Y to an assigned well f2j, if the approximate location 
of its minimum mi is known. 

Indeed, the function a{x) given, the global dynamics of Y is determined by the values a(mj), 1 < i < n, 
and the cooling rate 6. Moreover, for our analysis we can freely choose both a{x) and 9. 
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Let a(x) be smooth and attain its unique global maximum at to.;. Then we have 

a(rrii) > maxa(mj). (5-1) 

For instance, one can take a(x) = (A — a)/(l + (x — mi) 2 ) + a for some < a < A < 2. Then we can choose 
9 > 0, such that 

a(rrii)6 > 1, whereas a(m,j)6 < 1 for j^i. (5-2) 

With this choice of parameters, as i — » +oo, the particle leaves any well Qj, j =/= i, in a finite time according 
to (|4.3|) . Moreover, since all transition probabilities in ()4.4j) are strictly positive, the probability to enter the 
well fii after a finite number of transitions between the wells Oj, j ^ i, equals 1. Finally, upon entering f2j, 
the particle gets trapped there due to (|4.5[) . 



5.2 Looking for the global minimum 

Let M be the (unique) unknown global minimum of the potential U. To make a Levy particle settle near it, 
we have to determine the appropriate variable stability index a(x) and the cooling rate 9 such that 

a(M)9 > 1, whereas a(m l )9 < 1, m, ^ M. (5.3) 

Let <p(w) be an arbitrary smooth monotone decreasing function on K, < a < tp(u) < A < 2. Then we set 

a(y) = <p(U(y)). (5.4) 

It is clear that 

a(M) > a{rrii) for all m, ^ M, (5.5) 

and we choose 9 to satisfy relations (|5.3p . The trapping of the particle near the global minimum M follows 
from the argument of the preceding section. 



5.3 The local minimum with maximal energy 

Analogously, we can determine the coordinate of the (unique) 'highest' local minimum to, i.e. such that 

U(m) = max C7(mj). (5.6) 

l<i<n 

In this case, we should take a(y) = ip(U(y)) with an arbitrary smooth monotone increasing function ijj, 
< a < ip{y) < A < 2, which leads to the inequality 

a(m) > a(w,i) for all rrii ^ m, (5-7) 

and the arguments of the previous sections justify the success of the search. 



5.4 A local minimum with certain energy 

Finally, we can perform a search not for the global minimum of U but for a local minimum wie satisfying 
the condition U{w,e) < E for some energy level E. In such a case we take a piece- wise constant stability 
index 

a(y) = \ A ' U ®* E > (5.8) 
y; [a, U{y)>E 1 0<a<A<2, K ' 

and a cooling rate 9 satisfying conditions A9 > 1 and a£ < 1. Consequently, if for some i, the minimum 
U(nii) lies below the threshold E, then near this minimum Y behaves like a jump-diffusion Z^ A ' driven 
by a Levy flights process L^ A \ and thus the Levy particle gets trapped due to (|4.5p . If the well minimum 
lies above the level E, a Levy particle behaves like a process Z^> and leaves this well in finite time due to 
relations ([4~3]) and ([4~4]) . 

If there are several wells with minima below E, the Levy particle settles down near one of them. 
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Running the search for decreasing energies E\ > E% > we can also estimate the ground energy level 
E* = U(M) and to determine the global minimum M. Analogously, one can determine a local minimum 
rriE satisfying conditions IHttie) > E or E\ < U(mE) < -^2- 

We emphasise that the search algorithm described in this section requires no particular information about 
the potential U: to determine the local minimum jtie we use three numbers a, A and 9, the energy level E, 
and values of U'(x), 



6 Numerical example 

In the preceding sections we gave a theoretic justification of the search algorithm in dimension one. It is 
clear that a d-dimensional case, d > 1, does not differ much. It suffices to consider isotropic (spherically 
symmetric) d-dimensional Levy flights with a jumping measure v x (dy) = \y\~ d ~ a ^dy, and to calculate new 
values of transition probabilities in l]2.8[l writing d instead of 1 in the integrals. 
Similarly to (|3.3|) . we generate a d-dimensional Markov chain 



yfe = yfc-i - V[/(y fe _i)/i + 



#(y*-0 



Kk< 



(6.1) 



with some initial value yo- Isotropic a-stable random vectors ££(y) are obtained as marginals of a sub- 
ordinated standard d-dimcnsional Brownian motion (e.g., see Examples 24.12 and 30.6 in |Sat99| ). More 
precisely, let S be an a/2-stable strictly positive random variable with the Fourier transform 



Eexp(iujS) = cxp [-\u>\ a/2 (l - i sgn(w) tan( — )) ) , uel, 



and the Laplace transform 



Ecxp(— uS) = exp 



.a/8 



u > 0. 



cos(^), 

Let W be a standard Gaussian vector independent of S with a characteristic function 

Eexp(i(w,W)) =exp(-^-\ uj£R d . 



(6.2) 
(6.3) 

(6.4) 



Then one obtains a standard isotropic d-dimcnsional a-stable random vector L^"* 1 with the Fourier transform 



Eexp fi(w,L (Q) )) = cxp 



exp 



R d \{0} 

K d ' 2 \T{- 



e *(«^-l-i( W> y)I D (y) 



\d+a 



(6.5) 



2« d+a j 



oj G 



as 



L (Q) = c(a,<t)VS- W, c(a,d) = — 

v 2 



TT^COS 



/Trax [r(-f)| 

V 4 / r(^a) 



l/o 



Finally the random increments in (|6.1[) can be calculated as 

££( y ) = / l V«(y)L( Q (y)). 



(6.6) 



(6.7) 



6.1 Five- well potential in M 2 

To illustrate the efficiency of the method, we consider a two-dimensional potential function U given by the 
formula 



U(y) = 



l - 



l 



1 



1.5 



1 + 0.05(^ + ^2-10)2) l + 0.05(( yi - 10) 2 + y!) l + 0.03((jy 1 + 10)2 + yl) 



1 + 0.05((yi - 5)2 + (2/ 2 + lO) 2 ) 1 + 0.1(( yi + 5) 2 + (j/ 2 + 10) 2 ) 



(6.8) 



(l + 0.0001(2/^ + 2/ 2 2 ) 12 ). 
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Figure 1: A typical random search path (yk)o<k<n of a Levy particle (1.) and the values of the potential 
function U(yk) (r.). Thick lines on the left figure denote the boundaries of the attraction domains fii of the 
stable points m^, i = 1, . . . , 5. 



The function U has five local minima m;, i = 1, . . . , 5, with the following coordinates and energy values: 



mi p 




-9.73, 


-0.11), 


E/(mi) s 


^ -0.85, 


m 2 p 




-0.09, 


9.64), 


f/(m 2 ) r 


i -0.44, 


m.'i R 


*( 


9.59, 


-0.37), 


C/(m 3 ) R 


i -0.54, 


m i p 


*( 


4.92, 


-9.89), 


f/(m 4 ) R 


i -1.46, 


m 5 r 


*(- 


-4.79, 


-9.79), 


C/(m 5 ) r 


i -0.79. 



(6.9) 



We perform a search for a local minimum m^ having the energy less than E = —1, i.e. such that U{vo.e) < 
— 1. In our example, = m 4 . According to Section we choose a piece- wise constant stability index 

a(y) = { 1 - 8 ' if ^ < - 1 ' (6.10) 
yy> \l.l, ifC/(y)>-l, 1 j 

and 6> = 0.75. 

We perform 100 simulations of the Markov chain (|6.ip for n = 2 • 10 6 , initial conditions yo distributed 
uniformly in the square [— 20, 20] 2 , the initial temperature A = 10 4 , and the time step h = 0.1. The global 
minimum 1x14 was determined in 96 cases (U(y n ) = —1.46). The local minimum 1113 was located twice, and 
each of the minima mi and m 2 once. A typical random path (yfc)o<fe<« on the plain and the corresponding 
values of the energy function U (yt ) are shown on Figure [TJ 

We emphasise, that in our simulations we used only values of the potential U without any additional 
information about its geometry 



7 Conclusion and discussion 

In this paper we presented a new stochastic algorithm for global optimisation. It allows to determine a global 
minimum of an unknown potential U with help of simulated annealing of non-Gaussian jump-diffusions driven 
by the so-called stable-like processes, or Levy flights with a variable stability index a(x). We have shown 
that choosing ot(x) in an appropriate way, we can force the Levy particle to settle in a neighbourhood of 
the global maximum of U . We note, that the non-constant behaviour of the stability index is crucial for 
the success of the search, and a similar algorithm with usual spatially homogeneous Levy flights, i.e. when 
a(x) — a>o, leads to quite different results, see }Pav06a| . 

Our method has the following advantages in comparison with the Gaussian simulated annealing considered 
in the introduction. First, the search of the global minimum is non-local, i.e. when the annealed process 
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leaves a potential well, it does not necessarily pass to one of the neighbouring wells, but with strictly positive 
probability can jump to any well. Moreover, the probability to jump into the deepest well is maximal, if 
this well is also spatially the biggest, which is observed in typical potential landscapes, see |Sch97j . We do 
not expect that our algorithm would effectively detect the so-called 'golf-hole' wells. Mean transition times 
between the wells increase as a power of the large parameter A or, cquivalently, the current time t. We can 
easily obtain theoretic estimates for a number of transitions between the wells before settling in the deepest 
well. These estimates follow from the analysis of a discrete time Markov chain with transition probabilities 

Pi j = q% (mi)) i<ll a(mi)) ]- 1 ,Pii = 0- 

Second, we have more freedom to choose the parameters of the system. Indeed, if the values of U(rrii) 
are not known, there is no method which helps to determine the cooling rate 9. (One has the same problem 
to determine 9 in (|1.2p in Gaussian case.) However, in our algorithm, 9 is chosen together with a variable 
stability index a{x). 

Third, our method allows to drive the Levy particle into any well whose location is approximately known. 
We can also determine a local minimum with energy below the certain given value. The choice of parameters 
in these regimes is independent of the geometry of the potential. Such search regimes arc not possible in the 
classical setting. 

Finally, the temperature decreases polynomially fast in time, i.e. ~ t~ e , and not logarithmic. This 
significantly increases the accuracy of empirical estimates for the local minima locations. 

Although the theoretic basis for the success of the search is established, many questions are still open. 
For example, we have to understand how to choose the optimal pair a{x) and 9, which minimises the search 
time. Indeed, if m, is not a global minimum, we can reduce the life time of the particle in the neighbourhood 
of rrii making a(mj) small. On the other hand, in this case, the process Y will tend to make very big jumps, 
and thus jump out to one of the peripheral wells. As a consequence, the search would be slow, if the global 
minimum of U is attained in one of the inner wells. Thus, the value of a(m,) should not be very small to 
exclude very big jumps, and should be well separated from a{M) to block trapping in the false well. The 
problem of very big jumps can be also avoided by consideration of truncated Levy flights with maximal jump 
size not exceeding the size of the search domain. However, in this case the simulation of random input can 
be more complicated. We shall address these and other questions in our further research. 

Acknowledgements. This work was supported by the DFG research project Stochastic Dynamics of 
Climate States. The author thanks P. Imkeller for stimulating discussions. 
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